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ABSTRACT 

We present a smoothed particle hydrodynamic (SPH) simulation that reproduces a galaxy that is a moderate 
facsimile of those observed. The primary failing point of previous simulations of disk formation, namely excessive 
transport of angular momentum from gas to dark matter, is ameliorated by the inclusion of a supernova feedback 
algorithm that allows energy to persist in the model ISM for a period corresponding to the lifetime of stellar 
associations. The inclusion of feedback leads to a disk at a redshift z = 0.52, with a specific angular momentum 
content within 10% of the value required to fit observations. An exponential fit to the disk baryon surface density 
gives a scale length within 17% of the theoretical value. Runs without feedback, with or without star formation, 
exhibit the drastic angular momentum transport observed elsewhere. 

Subject headings: galaxies: formation, hydrodynamics, methods: N-body simulations 
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1. INTRODUCTION 

A fundamental assumption of the analytic work on disk for- 
mation of Fall & Efstathiou (1980) is that during protogalatic 
collapse the baryonic component conserves its specific angular 
mome ntum (AM). Pioneerin g N-body hydrodynamic simula- 
tions (Navarro & Benz 1991) showed that baryon cores, which 



develop naturally within CDM simulations, systematically lose 
specific AM to dark matter halos during merger events, con- 
trary to the assumption of Fall & Efstathiou. As detailed in 
White (1994), the resolution of this perceived 'numerical prob- 
lem' is widely believed to be the inclusion of feedback from su- 
pernovae (although for an alternative view see M ac Low & Fer- 
rara 1999), which is ne cessary in CDM cosm ologies to avoid 
the cooling catastrophe ( |White & Frenk 1991 ). To date, simu- 
lations incorporating feedback — using a variety of algorithms to 
model the process — have been largely u nsuccessful in preserv- 
ing the spec ific AM of the baryons (e.g. |Katz 1992|, N avarro & 
White 1993, [Navarro & Steinmetz 2000| ). Models, in which 
the gas is a rtificially prevented from cool ing before a 'cool- 
ing epoch' (Weil, Eke & Efstathiou 1998) have been consid- 



erably more successful. Alternatively, Dominguez-Tenreiro et 
al. (1998), argue that the inclusion of star formation can help 
stabilize disks against bar formation, and subsequent AM loss. 
Sommer-Larsen & Dolgov (2000) have shown that warm dark 
matter (WDM) does not suffer as significant an AM deficit 
problem due to the reduction in short-scale power compared 
with CDM. 

This Letter presents results of simulations which model star 
formation and feedback using an algorithm tested in detail 
in Thacker & Couchman (2000 — TC00). We achieve signif- 
icant success in reproducing fundamental properties of ob- 
served galaxies: the model galaxy does not suffer a catastrophic 
(^80%) loss of AM relative to the dark matter and has a disk 
scale length comparable with those observed and in accord w ith 
theoretically predicted values (e.g. Mo, Mao & White 1998). 



tion. The result of limited numerical resolution is a large dis- 
parity between the minimum simulation timescale and the true 
physical timescale associated with supernova feedback. Physi- 
cally, following a supernova, the gas cooling time very rapidly 
increases to (9(10 7 ) years over the time it takes the shock front 
to propagate. Thus it is the timescale — t^ — for the cooling time 
to increase markedly which is of critical importance. In the 
SPH simulation, "stars" form and "supernova feedback" occurs 
in gas cores which are typically less than h m j„/5 in diameter 
(h m i„ being the effective spatial resolution). Any rearrangement 
of the particles in such a small region leads to a density es- 
timate varying by at most ~ 7%. Since the cooling time in 
the model is dependent on the local SPH density, any energy 
deposited into a feedback region only reduces the local cool- 
ing rate by increasing the temperature; the SPH density cannot 
respond on a timescale which in real supernova events would 
very rapidly increase the cooling time. A successful feedback 
algorithm in an SPH model must thus overcome the fact that 
Psph does not change on the same time-scale as tg and conse- 
quently that cooling times for dense gas cores at T < 10 6 5 K 
and p > 2tiB cm" 3 remain short (< 10 Myr, as emphasized in 



Sommer-Larsen, Gelato & Vedel 1998) 



2. NUMERICAL LIMITATIONS OF FEEDBACK MODELS 



3. ALGORITHM AND INITIAL CONDITIONS 

We model a standard Cold Dark Matter model: £1/, = 0.1, 
flcDM = 0.9, h = 0.5, (78 = 0.6 and shape parameter T = 0.41 with 
the adiabatic Bond & Efstathiou (1984) CDM power spectrum. 
The same candidate halo as discussed in Thacker (1999) was 
selected for resimulation from a 100 3 dark-matter-only simula- 
tion of comoving width 48 Mpc. The halo, of mass 1 .66 x 10 12 
Mq lies on a filament ~ 2 Mpc long, and does not have a vio- 
lent merger history at low resolution. 

Three SPH simulations were run, one with no star formation 
or feedback (NSF), one with star formation but no feedback 
(NF) and one with both star formation and feedback (ES: En- 
ergy Smoothing, see below). The initial conditions for each 
simulation were prepared by creating four mass hierarchies in 
radial shells within the 48 3 Mpc 3 (comoving) simulation vol- 
ume. The central high resolution region is 6 Mpc in comoving 



It is not currently possible to model individual star formation 
and feedback events in cosmological simulations of disk forma- 
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FIG. 1 . — Density-temperature rilling functions for the NF run (left) and ES (right) at z=3 (top) and z=l (bottom). The magenta lines separate gas which can cool 
(below the line) and cannot cool (above the line) by z = at each epoch (t coo i = tg —t). 



diameter and has a total mass of 7.8 x 10 12 M Q and a parti- 
cle number of 2x65,454 yielding particle masses of 1.2 x 10 7 
M and 1.1 x 10 8 M Q for gas and dark matter respectively. 
The minimum 'glob' and dark halo mass resolutions are 52 
times higher, corresponding to the number of SPH neighbors. 
The dark matter particle mass is small enough to avoid spuri- 
ous t wo-body heating of the gas particles ( Steinmetz & White 
1997). Only the high resolution region includes SPH particles, 
which were given an initial temperature of 100 K. The particle 
positions were drawn from a 'glass' and the power spectrum 
was truncated at the Nyquist frequency of each hierarchy. A 
Plummer softening length of e = 3 kpc was chosen and the min- 
imum SPH smoothing length was /i m ,„ = 3.52 kpc. 

The star formation and feedback algorithms are described in 
detail in TC00, the key features are summarized here. Radiative 
cooling is included for a gas of metallicity 0.01 Z Q . The star 
formati on rate i s calculated using a Lagrangian Schmidt Law: 
M* = yjA-nGpg c*M g = c*M g /tf ree f a u. The SFR normalization 
(or star formation efficiency) was set to c* = 0.06, 50% higher 
than the local group estimate of Gnedin (2000). Each gas par- 
ticle carries an associated 'star mass' and can spawn two star 
particles each of half the mass of the original gas particle. This 
procedure leads to a delay between star formation and the asso- 
ciated feedback (see TC00 for a full discussion). After spawn- 
ing the first star particle, the gas particle mass is decremented 
accordingly. The Energy Smoothing algorithm of TC00 (ESa) 
smoothes 5 x 10 15 erg g" 1 of feedback energy over the SPH 
neighbor particles after a star formation event and allows this 
energy to persist in an adiabatic state for a time t\n = 30 Myr, 
after which the energy is allowed to radiate away. The 30 Myr 
period is motivated to coin cide with the lifetime of stellar as- 
sociations ( |Gerritsen 1997 ), and is longer than the 5 Myr value 
used in TC00. Star formation is only allowed in regions where 
the baryons are partially self gravitating, p\, > 0.2 pdm, and the 
local flow is converging (V.v < 0). 

It is worth noting that since the adiabatic period in the ES 



algorithm, t\ n, is not a function of resolution, the (resolution 
dependent) Courant condition does not guarantee a sufficiently 
short time-step for feedback regions to expand. We have found 
that dt < O.lti/2 is necessary. This criterion imposes a maxi- 
mum time-step length of 2 Myr although in practice the time- 
step is limited by other criteria. At lower resolution this cri- 
terion would become significant, which may partially explain 
the comparatively poor results for the ES feedback algorithm 
observed in TC00, where the algorithm was observed to have 
only a minor impact on AM values. 



4. RESULTS 

The simulations were integrated to z = 0.52, which re- 
quired 15,930(15,900;17440) timesteps for the ES(NF;NSF) 
runs respectively, requiring between 9 and 14 days on a 4 
node AlphaServer ES40. Minimum time-step values were 
dt m in =0.12 Myr for all simulations, and average values were 
dt =0.44(0.44;0.40) Myr. At z = 0.52, the virial radius was 
r2oo = 203.4 ± 0.4 kpc, and particles from the second mass hi- 
erarchy had reached 1 .5 r2oo preventing reliable further integra- 
tion. The relaxation argument in Thomas & Couchman (1992) 
implies that 4 particles per softening volume is necessary to 
avoid two-body heating in the dark matter at z = 0.52, which 
is satisfied in the densest regions of the simulation. At this 
epoch, there were 14,546(14,466; 14,403), 8,350(5,805;12,713) 
and 13, 193(17,264;0) dark matter, gas and star particles within 
r2oo- At this SPH resolution shock processes should be repre- 
sented with moderate accuracy within the halo ( Steinmetz & 
Muller 1993). Halo, disk and bulge masses are summarized in 
Table 1. Disk edges are found using a surface density cut as 
in TC00, with a slightly lower limit of 7 x 10 12 M Q Mpc" 2 due 
to the increased mass resolution. At z = 0.52, the disk edges 
defined in this way were at 29(19;13)kpc. A comparison of the 
gas content shows that at this epoch the ES galaxy has 3.1 times 
more gas than the NF run. 



Thacker & Couchman 



3 



Run 


M h 


M b 


M d 


B:D 




(10 12 M ) 


(10 1() M ) 


(10 10 M ) 




ES 


1.56 


6.69 


3.94 


1.69:1 


NF 


1.55 


6.54 


3.26 


2.01:1 


NSF 


1.55 


7.98 


2.15 


3.71:1 



TABLE 1. — Halo, disk and bulge masses, and bulge:disk ratios for the 
simulations. 

Although the first gas cores are in place by z = 9, star for- 
mation begins at z = 5.6 and the first feedback events occur at 
z = 3.95. A peak SFR of 51 M Q yr _1 is reached for the NF 
run at z = 2.75, while for the ES run the peak is 27 M Q yr" 1 at 
z = 3.42. Due to feedback the ES SFR is suppressed by up to 35 
Mq yr" 1 compared to the NF run between z = 3.4 and z = 1 .6. 
At z = 0.52 the average SFRs over the previous 0.4 Gyr are 1 .67 
M Q yr" 1 for the ES run and 0.96 M yr" 1 for the NF run. By 
z = 1 in the ES run 98% of the bulge mass is in place and for- 
mation of the extended gas disk begins largely after this epoch. 
The morphology of the ES galaxy is type SO, i.e. the system 
has a distinct disk and bulge component while the NF galaxy is 
closer to E7, with a thin gaseous disk embedded within it. 
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FIG. 2. — Specific AM for various components of the simulations 
at z=0.52. The ES simulation achieves a j value close to the com- 
bined observational and theoretical constraint jd ~ 0.82 ji,. Boxes are 
inferred from data assembled by Fall (1983). Scale lengths are de- 
rived from the exponential fits to the disks described in the text and 
Cj = jd I ji, values are also given. 

4. 1 . Effect of Feedback on Overcooling 

In figure[j] we plot the mass filling factor for the gas f m (n, T), 
defined as the fraction of mass per unit logarithmic interval at 
density n and temperature T. The integrated, constant density, 
cooling curves are also shown to separate gas which can cool 
from that which cannot. The z=3 plot for the NF run shows hot, 
low-density (halo) gas (T ~ 10 6 K, tig ~ 10~ 5 cm" 3 ), cooled gas 
cores (T ~ 10 4 K, rig > 10~ 3 cm" 3 ), and cold void gas (T < 10 4 
K, jib < 10~ 6 cm" 3 ). The ES simulation develops a new hot, 
high-density phase, which at z = 3 has T ~ 10 6 K, rag ~ 0.01 



(about the same as the hot halo) of the gas is in this phase. The 
'bridge' (/,„ > -2.0) between the hot-high-density phase and 
the hot-low-density phase (to the left of the cooling line) shows 
that some of the gas is being heated sufficiently and rising high 
enough in the halos to have t coo \ > 5 Gyr. This highlights the 
point that gas does not have to be blown-away to make it un- 
available for disk formation. However, it is clear from the z = 1 
plot that the predominant evolutionary pathway for feedback 
is recycling of gas back into the cold, high-density phase for 
which cooling times remain less than 1 Gyr. 

4.2. Angular momentum and disk scale length 

Navarro & Steinmetz (2000) argue, following work by Mo 
et al. (1998) and using the observations of Mathewson et 
al. (1992) and Courteau (1998), that the ratio of specific AM 
between disk and halo components, c,, should be given by, 
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at z = 0, where A is the average overdensity within the virial 
radius, V c is the disk circular speed, V200 is the circular speed at 
the virial radius, and A is the spin parameter of the halo. Thus, 
if a disk galaxy within a halo of overdensity 200 has the same 
V c as V200 and A = 0.05, then the disk component preserves one 
half of its initial specific AM. For the ES simulation, V200 =191 
km s" 1 , A = 0.075, and V c = 258 km s" 1 at r = 40 kpc which, 
extrapolated to z = 0, predicts Ct = 0.55. Parameters for the NF 
run are identical (to within 2%). 

In figure 2 we plot the specific specific AM of the halo, 
galaxy (bulge+disk) and dense cores at z = 0.52 for both sim- 
ulations. The most significant result is that the specific AM of 
the entire ES galaxy is only 10% lower than the predicted value 
of Cj = 0.55. Both the NF and NSF runs have minimal amounts 
of AM. This does not agree with the findings of Dominguez- 
Tenreiro et al. (1998), who argue that star formation can help 
prevent the formation of an angular-momentum-robbing bar. 
However, since the NSF simulation does not form a bar, our 
argument is not conclusive. As expected, the bulge compo- 
nents (not plotted) exhibit an extreme deficit of specific AM 
Ubulge ~ 0.01 j/,) having been formed primarily from first gen- 
eration infall that has been subjected to little feedback (even in 
the ES run). The dense cores, defined as all the baryons within 
r2oo satisfying 5b > 2000 (i.e. the galaxy plus satellites), have a 
high specific AM value because they include the contribution of 
one large core at 0.8r2oo (160 kpc) that heavily biases the r x v 
calculation. A large fraction of this AM will be shed as the core 
merges with the central galaxy. The presence of this satellite 
highlights the argument of Binney, Ortwin & Silk (2001): if 
one can expel the low angular momentum core, then baryons 
falling in from large radii, which naturally have a large specific 
AM, will produce a galaxy with a large specific AM value. 

Mo et al. (1998) have constructed an analytic model of disk 
formation that assumes a negligible exponential disk embedded 
within an isothermal dark halo. They parameterize the fraction 
of specific AM retained within the baryons in terms of cj and 
find a disk scale length given by, 



Rh = 



XGM, 



3/2 



2V 2 oo|£/ 1 | 1 / 2 



(2) 



For our simulated halos \Ei,\ = 3.45 x 10 J which implies 
cm" 3 , and log/,,, ~ -1.3 indicating that only a small amount R d = 13.2c, kpc. Using the cj values in figure 2, the ES disk 
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has an expected scale length of ^6.5 kpc, while for the NF 
and NSF runs the value is ~ 2.5 kpc. Least-squares exponen- 
tial fits to our disks (measured from e to the disk edge) give 
fl rf =7.6(4.7;3.6)kpc. 
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FIG. 3. — Tangential velocity and radial velocity dispersion for the 
ES galaxy at z=0.52, averaged within 104 particle Lagrangian bins. 
The tangential velocities compare well with the expected circular 
speed. 

4.3. Rotation curve and density profiles 

The large size of the stellar bulge makes an accurate deter- 
mination of the stellar disk rotation curve difficult. As a com- 
promise we concentrate on the gas disk rotation curve where 
star formation is ongoing. In figure 3 we show the tangential 
velocity of the gas in the ES run, compared to the expected ro- 
tation curve, which has been softened using the S2 softening 
shape (see TC00 for details). Rotation curves for the runs are 
broadly similar, peaking at 327(312;340)km s" 1 . The NF run 
has the lowest central density as the rapid creation of stars dur- 
ing the formation process results in a comparatively diffuse cen- 
tral core, which is an unavoidable artifact of the Tow' resolution 
of our models. Radial dark matter density profiles are almost 
identical, with p(r) oc r~ 2 interior to r = 6 kpc, similar to the re- 
sults of Tissera & Dominguez-Tenreiro (1998). A run without 
baryons has a profile that matches the Moore et al. (1999) fit. 

5. SUMMARY AND DISCUSSION 

Attempts to model star formation and feedback in a consis- 
tent way in previous cosmological simulations of galaxy forma- 
tion have had limited success because the minimum simulation 
timescale far exceeds the characteristic physical feedback time. 
This mismatch prevents feedback energy coupling effectively 



to the system; energy which is therefore unavailable to regulate 
star formation in the merging hierarchy. We have shown that 
by including a model for supernovae feedback in our simula- 
tion which allows energy to persist for a timescale of 30 Myr, 
we have overcome this obstacle. The technique has resulted 
in the preservation of half of the specific AM content of the 
baryons in the galaxy during collapse, 10% lower than the value 
required to reproduce observed disks at z = 0. (Some part of this 
deficit may be viewed as being due to the high Vc/Vuqq ratio in 
sCDM.) The disk scale length in the run with feedback is found 
to be within 17% of that predicted on theoretical grounds as- 
suming the same fraction of specific AM loss. The inclusion 
of star formation alone does not help the AM problem for this 
halo. Further, without feedback, there is a rapid exhaustion, by 
z ~ 0.5, of the gas supply for star formation: the overcooling 
problem in CDM cosmologies. Reducing the star formation ef- 
ficiency is not a solution as the specific AM values will remain 
similar to the NSF run. 

The effects of varying resolution are difficult to conjec- 
ture. The assumption of the appropriate — but fixed — physical 
time-scales is a barrier to the development of a resolution- 
independent model. It is worth noting, however, that our feed- 
back prescription will always produce feedback regions of a 
characteristic temperature between 10 5 and 10 6 K (see TC00). 
At higher resolution the first halos form with lower escape ve- 
locities and are consequently more susceptible to the effects of 
feedback, perhaps helping reduce the high central mass concen- 
tration evident in figure 3. Authoritative answers await higher 
resolution simulations. We plan to conduct a simulation with 
twice the linear resolution and eight times the mass resolu- 
tion. We are also conducting a 'survey' simulation using this 
feedback model at lower resolution, comparable to Evrard et 
al. (1994). 

The galaxy that we have simulated is large and represents 
a relatively rare galactic event and it is unwise to extrapolate 
the success of this particular model halo to a more general en- 
dorsement of CDM cosmogonies. The key result of this paper, 
however, is that hierarchical structure formation in not an anti- 
requisite for the successful formation of discs: the adoption of a 
plausible physical model for feedback can indeed regulate star 
formation and avoid catastrophic angular momentum loss. 
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